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We report results from high-resolution particle-mesh (PM) N-body simulations of 
structure formation in an O = 1 cosmological model with a mixture of Cold plus Hot 
Dark Matter (C+HDM) having Ocoid = 0.6, Qi, = 0.3, and Obaryon = 0.1. We present 
analytic fits to the C-fHDM power spectra for both cold and hot (z/) components, which 
provide initial conditions for our nonlinear simulations. In order to sample the neutrino 
O ' velocities adequately, these simulations included six times as many neutrino particles as 
O ' cold particles. Our simulation boxes were 14, 50, and 200 Mpc cubes (with Hq = 50 km 
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s ^ Mpc ^); we also did comparison simulations for Cold Dark Matter (CDM) in a 50 Mpc 
box. 

C-I-HDM with linear bias factor b = 1.5 is consistent both with the COBE data and 
with the galaxy correlations we calculate. We find the number of halos as a function of 



c§ ! mass and redshift in our simulations; our results for both CDM and C-fHDM are well fit 
by a Press-Schechter model. The number density of galaxy-mass halos is smaller than for 
^ . CDM, especially at redshift z > 2, but the numbers of cluster-mass halos are comparable. 
■ We also find that on galaxy scales the neutrino velocities and fiatter power spectrum in 
C-I-HDM result in galaxy pairwise velocities that are in good agreement with the data, and 
about 30% smaller than in CDM with the same biasing factor. On scales of several tens 
of Mpc, the C-f-HDM streaming velocities are considerably larger than CDM. As a result, 
the "cosmic Mach number" in C-I-HDM is about a factor of two larger than in CDM, and 
probably in better agreement with observations. 

Thus C-I-HDM looks promising as a model of structure formation. The presence 
of a hot component requires the introduction of a single additional parameter beyond 
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standard CDM — the light neutrino mass, or equivalently Qi, — and allows this model 
to fit essentially all the available cosmological data remarkably well. The tau neutrino is 
predicted to have mass of about 7 eV, compatible with the MSW explanation of the solar 
neutrino data together with a long-popular particle physics model. We outline a number 
of additional tests to which the C+HDM model should be subjected. 

Subject headings: cosmology: theory — dark matter — large scale structure of the universe 
— galaxies: formation — galaxies: clustering 

To be published in The Astrophysical Journal, October 1993. 
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1. Introduction 



The cold dark matter (CDM) model is perhaps the simplest potentially viable model 
for dark matter and structure formation in the universe. As is well known, CDM is based on 
the following set of assumptions (Blumenthal et al. 1984): the dark matter is cold, and the 
initial fluctuations are adiabatic, Gaussian, and have a Zel'dovich spectrum, as predicted 
by the simplest models of inflation. The "standard" CDM model makes the additional 
assumption that f2 = 1, with ilbaryon ~ 0.05 from standard big bang nucleosynthesis 
(Walker et al. 1991) and ficoid = ^ — ^^baryon; fitting observational data then requires that 
the Hubble parameter h = Hq/IOO km s~-^ Mpc ~^ 0.5 and that galaxy formation be 
"biased" (Davis et al. 1985). This attractive model had great success, but it is now well 
known to be in varying degrees of difficulty with many sets of observational data (see e.g. 
Davis et al. 1992 and references therein). 

Perhaps the simplest variant of CDM that remains viable has f2 0.2 with h ^ 1 
and a cosmological constant A = A/3Hq = 1 — Q ior consistency with inflation and with 
CMB constraints. This model is claimed to be consistent with the APM galaxy angular 
correlation function Wg{6) (Efstathiou, Sutherland, & Maddox 1990), with the observed 
rich cluster correlation function Cc{f) (Holtzman & Primack 1993), mass function (Lilje 
1992, Bacall & Cen 1992), and power spectra from clusters (Scaramella 1992), the CfA 
slices (Vogeley et al. 1992), and the Southern Sky redshift survey (Park, Gott, & da 
Costa 1992). However, there are several indications that Q Ri 1, for example CMB dipole 
vs. QDOT/IRAS data (Rowan- Robinson et al. 1990, Strauss et al. 1992), comparison of 
IRAS density and galaxy peculiar velocity data (Kaiser et al. 1991, Dekel et al. 1992), 
reconstructing Gaussian initial conditions from the POTENT analysis of galaxy peculiar 
velocity data (Nusser and Dekel 1992), and void outflow (Dekel & Rees 1992). While 
this evidence that f2 = 1 is still not compelling, and the arguments for a large Hubble 
parameter and an old universe do point toward smaller Q, we will adopt f2 = 1 and 
h 0.5 as hypotheses in this paper. 

The question arises whether any = 1 model with a physically motivated smooth 
spectrum of adiabatic Gaussian fluctuations can account for all the data now available, 
including the COBE CMB fluctuations (corresponding to scales of 3000-300 h~^ Mpc), 
large scale structure data (300-10 h"-*^ Mpc scales: galaxy angular correlations Wg{9), the 
cluster correlation function ^c, and galaxy streaming velocities), and smaller scale structure 
data (10 h~^ Mpc-10 h~-^ kpc: galaxy formation, correlations, and velocities)? 

One variant of standard CDM that has received much attention recently (e.g. Cen et 
al. 1992, Adams et al. 1993, Liddle & Lyth 1992b) keeps all the usual assumptions except 
the Zel'dovich primordial spectrum oc k"^ with n = 1, substituting instead "tilted" 
spectra with n ~ 0.5 — 0.7 which arise from more or less complicated inflationary models. 
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Such models have the virtue of being very well specified, with n being the only additional 
parameter beyond those of standard CDM. However, it appears that "tilted" CDM is 
marginal at best. For example, for n sufficiently small to account for the observed large 
scale structure, there is probably too little early galaxy formation. Of course, it is possible 
to get much more general non-Zel'dovich primordial fluctuation spectra from inflation (see 
e.g. Primack 1991 for a review), but these "designer spectra" are neither well motivated 
nor well specifled. 

Cold + Hot Dark Matter (C+HDM) is physically at least as well motivated as tilted 
CDM or any other variant of CDM that we know. Moreover it is well specified and has only 
one additional parameter beyond those of standard CDM, the neutrino mass or The 
currently available solar neutrino data suggest that the electron neutrino and one other 
neutrino, say the muon neutrino, have nonvanishing mass, with m(z^p) = (2 — 3) x 10~^ 
eV. In simple "seesaw" models of neutrino masses in which the three generations of heavy 
right-handed neutrinos are assumed to be degenerate in mass, m(z/T-) ~ 0.3{mt/mc)'^m{u^), 
where nit and nic are the top and charm quark masses, so m(z/r) should lie in the relevant 
mass range ~ 10 eV for C+HDM (for relevant references and recent examples of such mod- 
els, see Ellis, Lopez, & Nanopoulos 1992 and Dimopoulos, Hall, & Raby 1992). Moreover, 
the CHORUS and NOMAD Vfj^v-r oscillation experiments now underway at CERN could 
see a signal within a few years if these neutrino mass models are right and the mixings are 
large enough. 

While it is true that C-I-HDM does have an extra knob to adjust compared to CDM, 
it also has a physical feature that is suggested by the data: the neutrinos provide an 
unclustered dark matter component on small scales, which could help explain why dynam- 
ical estimates give O < 1 on small scales. The out-of-equilibrium relativistic Fermi-Dirac 
statistics of the neutrinos (once the neutrinos decouple, their momenta just redshift; see 
e.g. Weinberg 1972, p. 535) enhances this effect, as we will discuss below. 

The main objection to C+HDM in principle is the apparent unliklihood of having 
two different dark matter components each making comparable contributions to the mass 
density. Although one of the earliest C+HDM papers (Shafl & Stecker 1984) proposed a 
particle physics model to account for this, it must be admitted that we are unaware of any 
such model that is attractive. However, the entire particle physics Standard Model begs 
for further explanation, so it does not disturb us to contemplate one more feature that, if 
valid, would call for a more fundamental justiflcation. 

The question that wc will consider here is whether C+HDM can be valid: can it 
account for the astronomical data? Basic properties of mixed dark matter models were 
worked out some time ago (Fang, Li, & Xiang 1984, Valdarnini & Bonometto 1985, Achilli, 
Occhionero, & Scaramella 1985); and the fact that C+HDM is a promising model for large 
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scale structure was established by several previous linear calculations (Holtzman 1989, 
Schaefer, Shafi, & Stecker 1989, van Dalen & Schaefer 1992, Schaefer & Shafi 1992, Taylor & 
Rowan-Robinson 1992, and Holtzman & Primack 1993). A simplified nonlinear calculation 
in a 14 Mpc box has been done by Davis, Summers, & Schlegel (1992), with the initial 
neutrino fluctuations set equal to zero. Gelb, Gradwohl, & and Prieman (1993) includes 
a cdm simulation in a 200 Mpc box starting from a linear fluctuation spectrum inspired 
by C+HDM. In the present paper we present the flrst detailed nonlinear calculations 
for C+HDM of which we are aware, with proper initial conditions, sufficiently many hot 
particles to sample velocity space adequately, and a careful analysis of dark matter and 
galaxy correlations and velocities with comparisons to the available data. 

This paper is organized as follows: §2 describes our calculation of the cold and hot 
fiuctuation spectra and gives fitting functions, §3 explains the details of our numerical 
techniques and spectrum normalization, and §4 explains our biasing scheme and galaxy 
finding algorithms. The remaining sections present our results: §5 the abundance of halos 
from oTir numerical simulations compared to a Press-Schechter calculation, then the matter 
and galaxy correlation functions (§6) and velocities (§7), and in §8 the density distribu- 
tion function. In §9 we summarize our results and the status of C-I-HDM as a model of 
cosmological structure formation. It fares remarkably well! 

2. The spectrum of fluctuations 

The initial spectrum of fluctuations for our C-I-HDM model is computed using linear 
theory. The linear calculations trace the evolution of fluctuations in flve components: 
baryons, photons, cold particles, massive hot particles, and massless hot particles. The 
relative densities of baryons, cold particles, and massive hot particles are here assumed 
to be Oedm = 0.6, = 0.3, Qb = 0.1, using h = Hq/IOO km s'^ Mpc = 0.5 for the 
Hubble parameter. The energy density of the photons is determined from the observed 
mean temperature of the microwave background, taken here to be 2.7 K. The massless hot 
particles are identifled with two species of massless neutrinos, which determines the energy 
density in the massless hot particles as compared with the density of the radiation. The 
massive hot particle is taken to be a light neutrino with a velocity distribution given by 
relativistic Fermi-Dirac statistics. 

For the initial conditions for the calculations, we have assumed adiabatic perturba- 
tions with a scale-invariant spectrum as predicted by standard theories of inflation. The 
linear calculations are performed using the techniques of Holtzman (1989). The coupled 
evolution of perturbations is computed from very early times until the present. The angu- 
lar distribution of the radiation perturbations is expanded in Legendre polynomials, and 
sufficient resolution (up to ~ 600 orders at the time of recombination) is kept so that small 
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angle anisotropies can be calculated. The momentum dependence of the massive neutrino 
perturbations is determined by separately calculating the evolution of these perturbations 
at 15 different comoving momenta. For the massive neutrinos, only the lowest two angu- 
lar moments of the perturbations contribute gravitationally, so only these moments are 
computed, using techniques similar to those described by Bond and Szalay (1983). This 
computation, which involves the solution of integro-differential equations, becomes lengthy 
at late times. However, once the massive neutrinos become nonrelativistic, the high order 
angular perturbations become small, and their subsequent evolution can be computed with 
the full angular dependence, which is described by a coupled set of differential equations. 
The switchover is made at the epoch {1 + z) ~ 220, when the rest mass of the neutrino 
dominates the total energy even for the highest momentum neutrinos. Additional time is 
saved in computing the small scale perturbations by setting the neutrino perturbations to 
zero at intermediate times when they contribute less than 0.5 percent of the total gravita- 
tional perturbation; once the neutrinos become nonrelativistic and their fluctuations begin 
to grow, their gravitational influence is turned on again. 

A qualitative discussion of the linear fluctuation spectrum is given by Holtzman (1989). 
In Figure 1 , we present the computed spectrum for our cosmological model along with some 
analytic flts (dashed curves). Figure la presents the fluctuations in the cold particles, while 
Figure lb presents the fluctuations in the neutrinos. Note that even at the current epoch 
(top curves in Figures lab) the amplitude of the neutrino fluctuations is different from that 
of the cold particles; this arises because the neutrino perturbations are damped by free 
streaming and Landau damping until the massive neutrinos become nonrelativistic, and 
the neutrino perturbations have not yet had enough time to grow to match the fluctuations 
in the cold particles. The fluctuation spectrum of the baryons quickly grows to match the 
spectrum of the cold particles after recombination; by z=25, the baryon fluctuations are 
within a few percent of the cold particle fluctuations, and in our dissipationless calculations 
we can treat them with sufficient accuracy by including them with the cold fluctuations. 

The power spectrum of total density is P{k, z) = (0.7-\/-Pcoid + 0.3-\/^)^ • The analytic 
flts shown in Figure 1 are given by 

X [l + 1.2/c^/2 _ + 347(1 - Va/5)k^/'^ - 18(1 - 0.32a'^)k'^] ^ , (1) 

where a = {l+z)~^ is the expansion parameter, and k is the wavenumber measured in units 
of Mpc~^. These approximations are valid for k < 30Mpc~^ and z < 25. The accuracy 
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(maximum deviation) of the Pcoid approximation is 5%, while the neutrino spectrum is 
accurate within 25%. (At z = 15, when our simulations start, the maximum error in the 
neutrino spectrum is 20%, at A; 0.06 Mpc~^; the next-largest error then is about 10% at 
A; 0.3. Note that a 20% error in the neutrino spectrum leads to a 6% error in the power 
spectrum of total density, comparable to the maximum error in the cold spectrum.) At 
high k, the leading term in the spectrum for the "cold" component is the term with A;^/^, 
the last term in the denominator providing a correction. 

Wright et al. (1992) introduced the quantity "Excess Power" 

_ ap{8h-^ Mpc ) 
= ^■\,{25h-^ Mpc ) 

as a measure of the shape of the spectrum, and pointed out that EP ^ 1.30 ± 0.15 is 
required to fit the APM w{6). (Here crp(r) is as usual the rms mass fiuctuation in a top- 
hat sphere of radius r calculated using P(/c, z).) EP is defined so that EP = 1 for standard 
CDM. For our C+HDM model, EP = 1.37, consistent with the APM data. 



3. N-body Simulations 

3.1 Code 

Numerical simulations were done using standard Particle-Mesh (PM) code (Hockney 
&; Eastwood 1981). The equations we actually solved are given by Kates, Kotok, & Klypin 
(1991). Most of results are based on three simulations with 256^ grid points of resolution 
of gravitational forces. Three other simulations with lower resolution (128^ grid points) 
were mainly used for tests. Each simulation had seven sets of 128^ particles. One set of 
particles represented "cold" particles (WIMPs/axions plus baryons), while the remaining 
six sets were used to simulate "hot" particles (neutrinos). The particles had different 
relative masses: each "cold" particle had a relative mass 0.7 and each "hot" particle had 
the mass 0.3/6 = 0.05. 



3.2 Thermal velocities of "hot" particles 

The six sets were arranged in three pairs, particles of each pair having random "ther- 
mal" velocities of exactly equal magnitude but pointing in opposite directions. The direc- 
tions of these "thermal" velocities were random. The magnitudes of the velocities were 
drawn from relativistic Fermi-Dirac statistics: 

dn(v) oc — -, (3) 
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where vo{z) = (1 + z)ckT^/m^c'^, c 



being the velocity of hght, T,^ = 1. 



95 K, and 



rur^c^ = 91.5n^h'^ eV = 6M{n^/0.3){h/0.5feV. 



(4) 



Note that vo{0) = 1 .2{mjj/7eY)~^ km and that the rms velocity Vrmsl^) = 3.596t'o(-2)- 
Thus at z = 10, for example, Vj-^is = 290 km s~^, certainly high enough to suppress u 
clustering on galaxy scales. These relatively large velocities account for the falloff of A at 
high k in Figure la. 

The reason to have many particles to simulate the neutrino distribution is that the 
motion of "hot" particles due only to the thermal velocities must not generate any spurious 
fluctuations. The arrangement of hot particles in pairs having oppositely directed velocities 
locally preserves the center of mass, thus more closely simulating "thermal" velocities 
of "hot" particles. With the center of mass preserved, the power spectrum of spurious 
fluctuations due to the discreteness is restricted to have dependence on the wavenumber 
k, which significantly reduces propagation of errors to large scales, while the amplitude at 
small scales is small thanks to the large number of particles. 

Test runs have been done to estimate the effects of the discreteness. The test simula- 
tions had 32"^ particles in each set moving in a grid of 64''^ cells. The size of the box was 
chosen to be 14 Mpc. No initial fluctuations were imposed on the particles. Thus "cold" 
particles initially had no velocities and were placed in a regular cubic grid. "Hot" particles 
had thermal velocities, which due to the discreteness effects produced fluctuations of den- 
sity and displacements of "cold" particles. Two simulations were run: one with six sets of 
"hot" particles as described above, another (like that of Davis, Summers, & Schlegel 1992) 
with only one set of "hot" particles. The latter simulation did not preserve the local center 
of mass and fluctuations were larger. We started both simulations at redshift 1 -|- 2; = 20. 
At the flnal moment, corresponding to redshift zero, the level of density fluctuations was 
•\/^(l) = 0.16 for the simulation with six sets and ^/^(l) = 0.38 for the simulation with 
one set, where ^(1) is the correlation function at one cell spacing (0.22 Mpc). Thus the 
fluctuations in the six-set simulation were satisfactorily small — a factor y/6 smaller than 
in the one-set simulation. The correlation function at one cell spacing is sensitive to very 
short scales. Another statistic, which is sensitive to longer waves, is the kinetic energy. In 
the absence of fluctuations, the kinetic energy must decay as a~^, where a is the expansion 
parameter. Because of spurious fluctuations induced by the discreteness, the actual kinetic 
energy is larger than the initial kinetic energy rescaled to the final moment. The ratio of 
the actual to the rescaled kinetic energies at 2; = was 1.05 for the six-set simulation and 
1.43 for the one-set simulation. Visual inspection of dot plots of "cold" particles at the 
final moment confirmed our conclusion that the simulation with six sets of "hot" particles 
mimics the neutrino thermal velocities with much better accuracy. 
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3.3 Initial conditions and normalization 

We used our analytical approximations (Equations (1)) for the "cold" and "hot" spec- 
tra. Initial positions and velocities of particles were set using the Zel'dovich (1970) approxi- 
mation as was first described for the three-dimensional case by Klypin & Shandarin (1983). 
We did not differentiate the velocity potential numerically as Efstathiou et al. (1985) did, 
because it generates excessive numerical noise at short scales. Instead, the displacement 
vector was simulated directly. Phases of fluctuations were exactly the same for "hot" and 
"cold" particles. When generating velocities of "hot" particles, the "thermal" component, 
as described above, was added to the velocity produced by the Zel'dovich approximation. 

The amplitude of fluctuations was normalized so our realizations are drawn from an 
ensemble producing the quadrupole in the angular fluctuations of the cosmic microwave 
background at the lljiK level measured by COBE (Smoot et al. 1992, whose power- 
spectrum-normalized quadrupole actually includes only I > 3 angular components). To 
be more specific, we present expressions for the large scale angular fluctuations, density 
fluctuations, and velocity fluctuations. If Pcoid(^) is the power spectrum of density fluc- 
tuations of "cold" particles, then the corresponding angular fluctuations a^^, rms density 
fluctuations dp, and rms 3D peculiar velocity a"„, are 
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A f dx 

-^3- / -^Pco\d{k)Ji+iM, x = kRH, (5) 

Jo ^ 



2 /'"'max a /""-max 

^p,coid = ^ / dkk'^Pcoidik), (T^,coid = ^ / dkP^o\d{k), (6) 

^min ^ fcmin 

where Rh = 6 x 10"^/i~^ Mpc is the radius of the horizon, /cmin and /cmax are the limits 
imposed by the finite size of the computational box, and Jn is the Bessel function. The 
rms angular fiuctuations due to all multipole components with given / is then 

2_(2[+l) 2 2/ + 1 A 
^« - 4^ ""'^ " 7r2Z(/ + 1) ' 

where ^ is a normalization of the spectrum at small k: Pcoid = -^k. In order to normalize 
the spectrum of fluctuations to given rms of the quadrupole Q2 we set the normalization 
constant to 

A = ^-^QlRh ■ (8) 

Note that the integral for the quadrupole Q2 converges at very long waves and this is the 
reason why we can use Pcoid instead of the spectrum of total fluctuations: the hot and cold 
spectra are identical in the long wavelength limit. 

A number of authors (Liddle & Lyth 1992a, Krauss & White 1992, Davis et al. 1992, 
Adams et al. 1992, Lucchin, Mattarese, and MoUerach 1992, and Souradeep & Sahni 1992) 
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have recently reminded us that the CMB fluctuations observed by COBE could include 
some contribution from tensor modes (i.e., inflation- induced gravity waves, Starobinsky 
1979) as well as scalar modes (density fluctuations). However, this is a relatively small 
effect in many models of inflation, comparable to the uncertainty in the COBE quadrupole 
including the cosmic variance, and we will ignore it here. Note that a larger tensor /scalar 
ratio in the large-angle CMB fluctuations would imply a smaller amplitude for the scalar 
fluctuation spectrum, and therefore a larger linear bias factor for our C+HDM model. 

To normalize the numerical simulations, we estimate cr^ and cr„ expected in a box 
of the same size L as our computational box. To make these estimates we set the limits 
of integration to k^in — (27r/L)/\/2 and /Cmax = (27r/L)A^/2, where N is the number of 
particles along an axis (128 in this case) and the factor \/2 takes into account the fact that 
each harmonic in the simulation represents a small cubic element in the phase space. (The 
factor actually depends on the shape of the spectrum.) Then we make a realization that 
gives the same a„ as obtained from numerical integration of the spectrum. Usually the 
spectrum of density fluctuations is already accurate within 10-15%. Finally we estimate 
the power spectrum of density fluctuations in the numerical simulation and compare it with 
the predicted one. A correction to the amplitude is made if necessary. The agreement of 
the power spectrum of density fluctuations in the numerical simulation with the prediction 
of the linear theory is the flnal normalization criterion. 

3.4 Simulations 

Sizes of the computational boxes for the C+HDM simulations were 14 Mpc, 50 Mpc, 
and 200 Mpc (the Hubble constant is assumed tohe H = 50 km Mpc~^, i.e., h = 0.5). 
The simulations were started at redshift 15 and run to redshift zero with a constant step 
Aa = 0.01 in the expansion parameter a. Each simulation had 128^ "cold" particles 
and 6 x 128^ "hot" particles moving in a grid of 256^ cells. Thus, the smallest resolved 
comoving scale was 0.055 Mpc, 0.195 Mpc, and 0.781 Mpc for the three simulations, 
and correspondingly the mass of a "cold" particle was 6.29 x lO^M©, 2.86 x IO^Mq and 
1.83 X IO^^Mq. 

As explained in the previous subsection, the simulations were normalized to the COBE 
quadrupole Q2 = 17fiK. This normalization corresponds to a biasing parameter b = 1.5, 
where 1/6 is deflned as the rms fluctuations of mass inside a sphere of radius 8h~^ Mpc . 
Another way of expressing the normalization is to give the rms peculiar velocity at some 
scale as predicted by the linear theory. The total (i.e., with no smoothing) rms peculiar 
velocity of "cold" particles is ay = 750 km s~^, and the rms velocity of a sphere with radius 

comparison, CDM with the same biasing parameter 
b = 1.5 predicts ay = 660 km s~'^ and V50 = 190 km s~^. (We compare these linear Vji 
predictions with the data in §9, below.) 
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A model with the CDM initial spectrum was also simulated. We used the spectrum 
given by Bardeen et al. (1986). The spectrum was normalized to give the same biasing 
parameter 6 = 1.5 as for our C+HDM simulations. This normalization corresponds to 
the rms of the quadrupole Q2 = 8.5//K, which is not compatible with the COBE DMR 
results. The size of the computational box was 50 Mpc. The initial phases of fluctuations, 
initial moment, and time-step were exactly the same as for the C+HDM simulation with 
the same box size. 

Figure 2 presents the distribution of matter in a rectangular slice of 10 Mpc thickness 
in the 200 Mpc box simulation. The slice was chosen to pass through the highest density 
concentration in the simulation. The density concentration is a rich galaxy cluster in the 
upper right quadrant of the slice. Figure 2a presents only a small (8%) random fraction 
of "cold" particles. The distribution of "cold" points for which the density is higher than 
Ptotai > 10(/Ototai) is showu in Figure 2b (1/6 of the points). The distribution of 1015 
"galaxies" in the slice is shown in Figure 2c. "Galaxies" were first preselected as maxima 
of density, and then those with mass larger than 10^^ Mq are plotted. The area of the circle 
presenting a "galaxy" is proportional to its mass. The distributions of "hot" and "cold" 
particles are significantly different at small scales even at the present time. Figure 3 shows 
the distribution of both components in a 220 kpc slice in the 14 Mpc simulation. Note 
that while all large condensations in the plot are seen both in "cold" and "hot" particles 
(like the large group in the left top corner, a few big "galaxies" in the filament, and the 
filament itself), small halos of "cold" particles do not find counterparts in "hot" particles 
(for example, the two "galaxies" in the void at the middle right). 

4. Biasing schemes and galaxy finding algorithms 

In general, the distribution of galaxies ( "luminous" matter) does not follow that of the 
dark matter. Thus we say that the galactic distribution is "biased". There are different 
notions hidden behind the term "bias". The simplest one is the linear biasing parameter 
6, which is just the normalization of the spectrum of fluctuations. For the C+HDM 
model normalized to the quadrupole 17 jiK the biasing parameter is b = 1.5, meaning 
that in a sphere of radius 8h~^ Mpc the level of density fluctuations as estimated by the 
linear theory is 1.5 times smaller than the level of fluctuations of the number of galaxies. 
Thus with this large scale normalization, the C+HDM model must be "biased." Another 
important notion is "physical bias" : luminosity density is not proportional to the total local 
density because the efficiency of galaxy and star formation can depend very nonlinearly 
on the density and can also depend on other aspects of the environment. Note that this 
could happen even for a model with linear biasing parameter 6 = 1. Physical bias is 
perhaps very complicated, involving many poorly understood aspects of galaxy formation. 
Although CDM can perhaps be rendered more consistent with the observed large scale 
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galaxy distribution and other data by invoking a physical bias mechanism (e.g., Babul 
& White 1991), this weakens the predictive power of the CDM model which was in our 
opinion one of its most attractive features. 

There are some assumptions that can simplify the situation. Thus if some particle 
never resided in a high density region, it cannot be a part of a galaxy. On the other hand, 
if a dense clump has the mass and the radius of a large galaxy, it probably is a galaxy. This 
is not necessarily true for an object the size of a dwarf galaxy, however, because supernova 
explosions could have swept away most of its gas, thus stopping star formation. 

Various algorithms have been used to find "galaxies" in numerical simulations. Cluster- 
analysis or "friends-of-friends" (Efstathiou et al. 1985, Prenk et al. 1988) has the advantage 
that no specific shape of an object is assumed. Unfortunately it also has significant disad- 
vantages. If the neighborhood radius (effectively a density threshold) is small, it does not 
find all halos (especially in low density regions). If the radius is larger, then it does not 
resolve halos in dense regions: in the case of an object which looks like a cluster of galaxies 
(judging from its radius and mass) the algorithm assigns all the mass of the cluster to one 
"galaxy", thus mixing together galaxies, groups, and clusters. 

We decided to use simpler prescriptions. Three different procedures were used to 
identify the "luminous" matter in our simulations. All these procedures deal with the 
density at a given moment of time. When possible, we compare results obtained with 
different methods or with different parameters in the same method. 

i) Density thresholding. If we do not need to deal with separate objects (for example, 
in the case of the correlation function), we find the density of "cold" particles on our 
standard 256^ mesh and consider a cell "luminous" if the density is larger than some 
threshold pthr,coid- The luminosity density is assumed to be directly proportional to the 
density at the cell. This simple prescription removes essentially all points from low density 
regions that have not yet collapsed. The density threshold is the free parameter. For the 
50 Mpc simulation the cold density threshold was chosen to be Pthr,coid = 25(pcoid)) which 
corresponds to the mass in a cell larger than 8.9 x 10^ Mq. This threshold is a reasonable 
compromise between expectations for the galactic mass, overdensity for collapsed objects, 
and statistical fluctuations in the correlation function. We discuss below the dependence of 
the correlation function on the threshold. The situation is worse for the 200 Mpc simulation 
because each "cold" point has enough mass to represent a galaxy. Nevertheless, we chose 
a threshold pthr,coid = 5(pcoid), thus selecting cells that collapsed in the simulation. This 
corresponds to the mass in a cell larger than 1.1 x IO^^Mq. To be precise, the thresholding 
of density does not guarantee that the object had undergone collapse. But in the case of the 
C-I-HDM model the effective slope of the power spectrum on galactic scales is about —2.5, 
which implies that the formation of an object goes first through the formation of pancakes 
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(e.g., Klypin & Melott 1992). The Zel'dovich approximation for pancake formation shows 
that the points with density at the threshold 5(pcoid) will collapse very soon. 

ii) Maxima of density. If we need to deal with separate objects (mass function, pairwise 
velocities), we find maxima of density. The total ("cold" plus "hot") density was produced 
on the 256^ grid. Then all local maxima above some threshold pthr are found. The way we 
assign mass to a maximum ("galaxy") depends on the particular problem we deal with. 
For counts of objects at different redshifts the mass of a "galaxy" is the mass in the cube of 
3x3x3 cells centered on the maximum. Effectively this corresponds to the radius of objects 
being 1.5 cell units. For simulations with 14 Mpc and 50 Mpc box sizes, this corresponds 
to 82 kpc and 293 kpc. For all the simulations the threshold of the density was chosen 
to be Pthr = 50(ptotai)- The Press- Schecter (1974) approximation and visual investigation 
of dot plots were used as guides for the choice of the parameters: with smaller size of the 
template we significantly underestimate the mass of halos in a dense environment, while 
with lower density threshold we overcount objects as compared with the Press-Schecter 
approximation. 

iii) Maxima of density inside a sphere. For the analysis of peculiar velocities of "galaxies" 
we do not need to know the mass of a "galaxy" precisely, but we would like to know its 
position slightly more accurately. After locating positions of maxima as for method (ii), 
we place a sphere of radius r at the position of each maximum and find the center of mass 
of all "cold" particles inside the sphere. We then displace the center of the sphere to the 
center of mass and repeat the procedure. After two iterations the positions of the spheres 
essentially stop changing. The typical distance between the final and initial centers of a 
sphere is less than 1/4 of the cell size, and the increase in the number of points inside a 
sphere is 20-40%. One can expect that this procedure gives a slight improvement over 
biasing scheme (ii). The typical radius of a sphere was 0.5 of the cell size. 

5. Abundance of halos 

In order to estimate the number density of galaxies at various redshifts in our sim- 
ulations we found all local maxima of the total density above the overdensity 50 (§4, 
prescription (ii)). Figure 4 shows the number density A^(> M) of dark halos with mass 
larger than M. The three bottom full curves (mass limits M = 3 x 10^^ Mq and higher) 
are for the simulation with 50 Mpc box size. The three top full curves correspond to the 
number of objects in the simulation with 14 Mpc box size. The dot-dashed curves show 
results for the Press-Schecter approximation, which were estimated as follows. For the 
Gaussian filter with radius Rf and mass M = (27r)"^/^po-R/ the number density of halos 
with mass larger than M is 

N{> M, z) = r nim, z)dm =^ T exp ( ) ^ , (9) 
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where 



eiz)= / k^P{k,z)e-^''''f^'dk/ / k'^P{k,z)e 



(10) 



and P{k, z) — (0.7\/Pcoid + 0.3\/^)^ is as usual the power spectrum of the total density. 
The parameter S^. was considered as a free parameter. The best results for 14 Mpc and 
50 Mpc boxes and for the CDM model simulated in a 50 Mpc box were found for 5c = 1-60, 
which is very close to the value given by the top-hat model. 

As one could expect, because long waves are not present in a finite box and because 
high peaks are rare in the box, the number of massive objects at high redshifts in the nu- 
merical simulations was below predictions of the Press-Schecter approximation. Although 
(as we explained in §4) the parameters of our galaxy finding algorithm were chosen to 
have results consistent with the Press-Schecter approximation, the fact that we succeeded 
for such a wide interval both in mass and redshift should be considered as a remarkable 
success for the approximation (cf. Efstathiou & Rees 1988, Bond et al. 1992). 

Figure 5 shows the number density of halos per unit volume and per unit logarithmic 
mass interval at redshift zero. The two full curves present the number of halos in C+HDM 
simulations for 50 Mpc box size (right curve) and for the 14 Mpc simulation (left curve). 
Results for the CDM simulation are shown by the dash-dotted curve. The dashed curve 
presents an analytical approximation: 



The power law slopes in this approximation were taken from the Press-Schecter approxi- 
mation for a slope n = —2.5 in the power spectrum, which is appropriate for the C-I-HDM 
model. We did not estimate the amplitude and the characteristic mass scale using the 
P-S approximation, so Equation (11) should be considered as a fit to the simulation data. 

It is interesting to compare the numbers of "galaxies" in the C+HDM and in the 
CDM simulations. The number density of "galaxies" in the CDM simulations at different 
redshifts is shown in Figure 6 (triangles are for the numerical results and dashed curves 
are for the Press-Schecter approximation with the parameter dc = 1.60). Full curves 
are the results of the Press-Schecter approximation for the C+HDM model. At redshift 
z = the CDM model predicts about twice the number of dark halos for masses M = 
10^^ Mq — 10^^ Mq . For less massive objects both theories predict essentially the same 
number of galaxies. The only masses where the C+HDM beats the CDM model are those 
of rich galaxy clusters M > 10^^ Mq . At higher redshifts the C+HDM model predicts 
progressively smaller number densities of objects. 





(11) 
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All these curves are for CDM and C+HDM normalized with the same linear bias factor 
b = 1.5. As is well known (Davis et al. 1985), CDM with this low a bias factor predicts 
galaxy peculiar velocities that are too large on small scales; as we will discuss shortly, we 
confirm this in our CDM simulations. CDM with b = 2.5 gives galaxy number densities 
that are more like our 6 = 1.5 C+HDM model. This is illustrated by the dot-dashed curve 
on the Figure. 

We note that the C+HDM model predicts quite substantial numbers of massive galax- 
ies at high redshifts. For example, at redshift z = S each cube of a size 100 Mpc contains 
f=:i 180 galaxies with mass larger than lO-'^-'^ Mq . Unfortunately, it is difficult to make a 
conclusive comparison with observations because there is still no reliable estimate of the 
number of galaxies at high redshifts with any given minimum estimated mass or circular 
velocity. The number density of QSOs is another test for the theory (Efstathiou & Rees 
1988); this has recently been considered in a linear calculation as a possible constraint on 
C+HDM (Schaefer and Shafi 1993). At redshift z = 2 the number of QSOs brighter than 
L = 10 erg s~^ is Nq ^ 10"^ QSOs per cubic Mpc {h = 0.5). If we assume that galaxies 
with mass larger than 3 x 10^^ Mq are responsible for the QSOs, then in the C+HDM 
model the fraction of galaxies hosting a QSO is about 2.5 x 10~^ and each of the host 
galaxies deposited about one per cent of its mass into the central mass associated with 
the QSO. Thus, the C+HDM model probably is compatible with the observed number of 
QSOs and galaxies at redshifts up to 2; = 2 — 3. But the C+HDM model will have diffi- 
culties if significant numbers of massive galaxies or bright QSOs are found at 2; > (4 — 5). 
For example, at 2; = 5 the number of halos in the C+HDM model with mass larger than 
10^^ Mq is 10"^ Mpc Thus if the number of QSOs per cubic Mpc remains Nq 10 
up to 2; = 5 (these QSOs could then plausibly generate enough UV radiation to account for 
the Gunn- Peterson test and the proximity effect, see Madau 1992), then every 10"*^^ Mq 
dark C+HDM halo should have a QSO inside. 

Estimates of the number of very rare high redshift objects based on the Press-Schecter 
approximation should be treated with caution, however. It is not clear whether the ap- 
proximation provides a good fit for very rare objects. And even in the framework of the 
Press-Schecter approximation there is a possibility to tune parameters to give many more 
dark halos at high redshifts. The parameter 5c, here assumed to be 1.60, could be smaller 
at higher z giving more objects. For example, Efstathiou & Rees (1988) used 6c = 1.33. 
At z = 4.5 the C+HDM model with 4 = 1-33 predicts A^(> 3 x 10^^) = 2.5 x 10"^, which 
is probably enough to account for observed high redshift QSOs. 

6. Correlation functions 

The evolution of the correlation function of the total density in the simulation with 
50 Mpc box size is shown in Figure 7 . The curves in the plot are labeled by the expansion 
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parameter. The correlation function ^(r) evolves much as in the HDM model: ^(r) is well 
approximated by a power law, with the slope increasing with time. At the final moment 
the slope of the correlation function is about 7 = 1.8, which is the observed slope of the 
galaxy correlation function, but the correlation length is 2.2h~^ Mpc and is too small as 
compared to the correlation length of galaxies. The dashed curve in this figure shows the 
correlation function of the "cold" matter at the same epoch. It is close to the correlation 
function of the total density, which indicates that there is not much difference between the 
distribution of "hot" and "cold" particles. For comparison we also present the correlation 
function of the dark matter in the CDM simulation. It is higher than that of the C+HDM 
simulations and it does not behave as much like a power law. 

The correlation length of the total matter is even smaller for the 14 Mpc box: 0.82h~^ 
Mpc. It is slightly larger in the 200 Mpc box simulation (Figure 8 ). The length is 2.9h~^ 
Mpc for this case, but it is still too small. The small correlation length of the dark as 
well as the "cold" matter indicates that the C+HDM model needs a bias to be consistent 
with the observed correlation function of galaxies. Normalizing to the COBE quadrupole 
led to the biasing parameter 6 = 1.5. If we assume that the correlation length for galaxies 
is 5.5/i~^ Mpc and take 2.9h~^ Mpc as an estimate of the correlation length in the 
C+HDM model, than the biasing parameter as found from the numerical simulations is 
b = \l (5.5/2.9)^-^ = 1.8, which is reasonably close to the COBE normalization biasing 
parameter. The actual biasing parameter for simulated "luminous" matter will be slightly 
higher (6 = 1.9), but still in reasonable agreement with the COBE normalization. 

To estimate the correlation function of "luminous" matter in the model, we use biasing 
algorithm (i) of §4 (density thresholding). Figure 8 shows the correlation functions of 
"galaxies" for the 50 Mpc C+HDM simulation (dashed curves, mass in each selected cell is 
larger than 8.9 x lO^M©) and the 200 Mpc simulation (full curves, mass in a cell is larger 
than 1.1 X lO^^M©). For comparison the correlation functions of the "cold" matter arc also 
presented (two bottom curves). The dot-dashed line is the power law ^ = (5.5/i~^/r)^'^, 
which gives a good approximation to the "galaxy" correlation function from 200 kpc to 
20 Mpc. 

The correlation function of "galaxies" depends on the biasing prescription. We found 
as expected that the correlation function is increased if regions with higher density or if 
more massive clumps are selected. The influence of different effects on the correlation 
function is demonstrated in Figure 9 . To estimate the correlation functions presented 
in the plot, we have not used the usual FFT technique, but instead we counted pairs of 
objects as a function of the distance between them. This accounts for small variations 
compared to previous plots. The pairs of "galaxies" were weighted by the product of the 
masses of the objects. This is needed in order to compensate the "overmerging" effect 
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in numerical simulations: because of lack of resolution and because no energy dissipation 
takes place when "galaxies" form, substructures disappear when objects of smaller size 
merge to produce a larger object. 

The two full curves in Figure 9a show how the correlation function scales with the 
density threshold in the 50 Mpc box simulation. The lower curve is for the density threshold 
10 and the upper curve is for the density threshold 25. The dashed and dot-dashed curves 
show results for "galaxies" found as maxima of number of "cold" particles inside a sphere 
of 100 kpc radius in the 50 Mpc box simulation. The correlation function of "galaxies" 
more massive than 3 x 10"^° Mq (913 objects) is shown as the dashed curve. The dot- 
dashed curve is for the 87 "galaxies" with mass larger than 3 x lO"^-*^ Mq . On scales below 
~ lh~^ Mpc , the variations in the "galaxy" correlation function are about a factor of two, 
the slope and the amplitude being higher for objects that are denser or more massive. On 
larger scales all curves are close to each other and the correlation length essentially does not 
depend on particular details of the biasing prescription and is equal to (4 ± 0.5)h~^ Mpc . 
For comparison we show results for galaxies in the CfA catalog (Davis & Peebles 1983). 
Although at small scales there is a reasonable agreement, "galaxies" in the simulation are 
systematically less correlated than in observations. The lack of long waves in the numerical 
simulation mainly accounts for the difference. If the long waves are included (as in the 200 
Mpc simulation), the correlation length for "galaxies" rises to the level (5 — Q)h~^ Mpc . 

Figure 9b shows effects of peculiar motions on the galaxy correlation function. Full 
curves in the plot present correlation functions in redshift space (peculiar velocity along 
the line of sight is interpreted as the additional displacement) for 50 Mpc (bottom curve) 
and 200 Mpc (top curve) simulations. Corresponding correlation functions of 913 and 2831 
"galaxies" in real space are shown as dashed curves. "Galaxies" in the 200 Mpc simulation 
were found as maxima of the number of "cold" particles inside a sphere of 300 kpc. The 
peculiar motions have two effects. Correlations at small scales are reduced because large 
velocities inside virialized groups spread the groups along the line of sight. This affects 
scales up to ~ 2h~^ Mpc , which corresponds to peculiar velocities of ~ 200 km s~^. The 
large-scale motions increase correlations, but the effect is small in our case. The markers 
on the plot present the correlation function in redshift space for galaxies in the CfA slices 
(Vogeley et al. 1992). 

7. Velocities 

One of the greatest difficulties for the CDM model is getting the velocities right on 
small scales. For the "unbiased" (i.e., 6=1) CDM model the peculiar velocities were much 
too large (Davis et al. 1985) compared to observed velocities in pairs of galaxies (Davis 
& Peebles 1983). Although "velocity bias" (Carlberg et al. 1990, Couchman & Carlberg 
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1992) perhaps reduces the pairwise velocities, it apparently does not entirely eliminate the 
contradiction with observations. Moreover, the amount of velocity bias, if any, apparently 
depends on the simulation technique and galaxy-finding prescription. As we discuss at the 
end of this section, our simulations do not support the suggestion (Carlberg et al. 1990) 
that dynamical friction leads to velocity bias. 

For the "biased" CDM model of Efstathiou et al. (1985) with biasing parameter h = 
2.5, the magnitude of the velocities was about correct, but the dependence of the velocities 
on the distance between galaxies was wrong: while the observed pairwise velocities slowly 
rise with the projected distance a = (340 ± 40)(/irMpc)°'^^'''° °^ km s~^ (Davis & Peebles 
1983), relative velocities of pairs of "galaxies" in the model decreased with distance for 
r < 1 Mpc and only then started to rise. 

In the C+HDM model the total rms velocity of matter relative to the rest frame is cr„ = 
750 km (the thermal velocities are not included), which is larger than ay = 660 km s~^ 
in the CDM model normalized to the same biasing parameter b= 1.5. However, most of the 
velocity in the C+HDM model is due to the long waves. For example, the rms velocity of a 
sphere with radius 50h~^ Mpc , which is the characteristic scale of the Great Attractor, for 
the C+HDM model is V50 = 360 km s~^ compared to the CDM value V50 = 190 km s~^ 
Because of nonlinear effects peculiar velocities in the numerical models were higher, but 
not much higher, than the linear theory predictions. For the 50 Mpc box simulation the 
linear theory predicts the rms 3D velocity 326 km s~^, while in the simulation it was 
384.6 km s~^. The 200 Mpc box simulation had larger velocities: (j„ = 623 km s~^ was 
predicted by the linear theory and almost the same value (617 km s~^) was actually found 
in the numerical simulations. 

The distribution of relative peculiar velocities of pairs is described by several functions. 
The first moment (^21) is the mean peculiar velocity along the line joining the pair (a 
positive value means that the two points move apart from each other faster than the 
Hubble fiow). The second central moment of the distribution of pair velocities has two 
components: the velocity dispersion along the separation vector (v|) and the velocity 
dispersion perpendicular to the vector {v'j_). Figures 10 and 11 present these components 
of the pair-wise velocity distribution in the 50 Mpc box simulations for both C+HDM 
and CDM models. The top panels show (f2i) for "cold" particles (dashed curves) and for 
"galaxies" (full curves). The dot-dashed curves stand for the velocity equal to the Hubble 
velocity: a pair on the curve would have zero proper velocity, thus indicating a stationary 
pair in a statistical sense. The bottom panels present the rms values of v\\ (full curves) 
and v± (dashed curves) for "cold" particles (the two top curves) and for "galaxies" (the 
two bottom curves). Because of the lack of resolution at small scales, one expects that the 
velocities in the simulation are smaller than they would be in reality. But the resolution 
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(one cell size) in the simulation corresponds to Q7h ^ kpc, so scales above 200 kpc-300 kpc 
were well resolved. 

The "galaxies" were found using algorithm (iii) of §4 (density maxima on the 256^ 
grid, then maxima of numbers of "cold" points inside a sphere of 100 kpc radius). The 
total number of "galaxies" was 913 in the 50 Mpc C+HDM simulation, and 1758 for 
the CDM simulation. The "galaxies" were more massive than 3 x 10-^° Mq . For the 
C+HDM simulation the rms peculiar velocity of the "galaxies" relative to the rest frame 
was 340 km s"-*^, which is only slightly smaller than the velocity of all "cold" particles 
(384.6 km s"-*^). But the pairwise velocity of "galaxies" differs significantly (by a factor 1.8 
at lh~^ Mpc ) from that of the dark matter at distances smaller than 3h~^ Mpc , while 
they are almost the same at larger scales. For the CDM simulation the situation is rather 
different. There is a larger difference between the rms velocities of "galaxies" and the dark 
matter relative to the rest frame (416 km for galaxies and 582.5 km s~^ for the dark 
matter). The difference of pair- wise velocities also is larger at lh~^ Mpc (by a factor of 
two), but at larger distances the difference does not go away. Even at 7h~^ Mpc the ratio 
of velocities is about 1.3. 

While (f2i) for "cold" particles in the C+HDM model shows that evolution still goes 
on in this population (slight outflow at scales < 1 Mpc and inflow at (1-2) Mpc), the 
pairs of "galaxies" at scales less than ^ 2h~^ Mpc look like stationary objects with their 
peculiar velocities cancelling the Hubble expansion. The same is true for the CDM model. 

In order to obtain pair-wise projected velocities, we place two "observers" at opposite 
corners of the computational box and simulate line-of-sight velocities Vp of "galaxies" by 
adding the line-of-sight component of the peculiar velocity to the distance scaled by the 
Hubble constant. The projected distance between a pair of objects is estimated in the 
usual way as Vp = tan(6'i2/2) (Vpi + Vp2)/H, where H is the Hubble constant and 6*12 is 
the angle between the objects. The distribution of velocity differences Vpi — Vp2 is then 
studied as a function of Vp. Following Davis & Peebles (1983), we do not include pairs with 
AV > 1000 km s~^ in the analysis. The solid curves in Figure 12 represent the velocity 
dispersion of "galaxies" along line of sight as a function of projected distance for the 50 Mpc 
box simulations for both the C+HDM and the CDM models. The dashed curves in the 
plot are for all "cold" particles and the dot-dashed curves show the law 3A0{hr)^-^^ km s"-*^, 
which is a good fit to the observations. The C+HDM model produces projected pair- wise 
velocities that give a very nice fit — both in amplitude and in the trend with the distance 
— to the observed velocities. The CDM model with the same bias factor {b = 1.5) gives a 
systematically larger amplitude by a factor 1.2-1.3. 

The dependence of the velocity dispersion on different effects in demonstrated in 
Figure 13 . The top panel shows results for the 200 Mpc box simulation. The velocity 



19 



dispersion is higher in this case, but in general the trends are the same and the difference 
in the amphtude is not too large. The fitting law is 380(/j.r)° -^^ km s~^, which is still 
compatible with the results of Davis & Peebles (1983). We found no dependence on the 
mass of galaxies, but with the available statistics the range in mass was not large. The 
long-dashed curve in the bottom panel gives results for 481 "galaxies" with mass larger 
than 10^^ Mq . If instead of cutting velocity differences at 1000 km we consider only 
pairs with the relative distance smaller than 10h~^ Mpc , the dispersion rises by only 
10 per cent (short-dashed curve). The dispersion significantly grows if the cutoff in the 
velocity difference is increased to 1500 km s~^ (dot-dashed curve). The same effect was 
found in the CfA catalog (Davis & Peebles 1983). 

The difference in the velocities of the dark matter particles and "galaxies" is called the 
velocity bias. We find that the global velocity bias measured as the ratio of rms velocities 
of all points to the rms velocities of "galaxies" is quite small — 1.13 for the 50 Mpc box 
simulation and 1.05 for the 200 Mpc simulation. We also have shown that at scales larger 
than 5h~^ Mpc, the velocity distribution of "galaxies" approaches that of "cold" particles. 
The difference appears on small scales, where velocities of "galaxies" are about half those 
of "cold" particles. The difference becomes even larger if we select those "cold" particles 
that reside in high density regions. Particles with p > 25(pmean) have velocities that on 
average are a factor 1.2 higher than the average for all "cold" particles. 

One source of the bias is the internal random velocities inside "galaxies" . The veloci- 
ties are present in the analysis of "cold" particles, but they are removed from the velocities 
of "galaxies". In our simulations internal velocities account for about half of the effect. 
Carlberg et al. (1990) suggested that dynamical friction is the main cause of the velocity 
bias they found. Because the acceleration due to the friction is proportional to the mass 
of an object, more massive objects will tend to sink to the bottom of the gravitational well 
and move slower. Thus, one expects that if dynamical friction were important, more mas- 
sive "galaxies" would move slower than less massive ones. Figure 14 shows the distribution 
of the total velocity of galaxies relative to the rest frame in the simulation with 50 Mpc 
box (913 objects). The curve on the plot presents the averaged rms velocity. Although 
the circular velocities of "galaxies" range from 50 km s~^ to 500 km s~^, there is no indi- 
cation that the velocity decreases with mass. (Note that with our force resolution in this 
simulation of 0.2 Mpc and mass resolution of about 3 x 10^ Mq , overmerging (see e.g. 
Gelb 1992) is not a serious problem for our galaxy velocity calculations. We have made 
plots of the density distribution in some of our groups and clusters, and checked that our 
galaxy-finding algorithm finds many "galaxies" there. Recent CDM calculations that also 
find no evidence for velocity bias include dissipationless simulations by Ueda, Itoh, & Suto 
1993, and hydrodynamic simulations by Katz, Hernquist, & Weinberg 1992.) 
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8. Density distribution function 

The density distribution function is an indicator of the dynamical state of a system. In 
the hnear regime the distribution function is a gaussian function. But even for a relatively 
small amplitude of the fluctuations, the distribution function significantly deviates from 
the gaussian law at large densities because rare high fluctuations become nonlinear and 
produce a small number of high density objects that cannot be found in the gaussian fleld. 

There are two asymptotics for the evolution of the density distribution function. For a 
model dominated by caustics (like the HDM model) the distribution function was found by 
Kofman (1991). The distribution function -P(p), which we define as a probability to find a 
volume element with the density p, has a tail P{p) oc for the pancake approximation. 
Another interesting case is a hierarchy of A/"— body correlation functions (Balian & Schaeffer 
1988), which probably arises at an advanced nonlinear stage. If the correlation functions 
scale as ^n{tii • • • , tn) = A'''(-'^~-'^)^jv(Ari, . . . , Ar^r), then the probability to find N points 
in a volume is Pjv N~'^'^^, < a; < 1. It was shown that for large N the probability 
Pjv has the scaling form ~ h{N/Nc)/{^Nc), where f is the average correlation function 
in the volume and Nc = nl^^, where n is the mean number density. The function 
h{x) decreases exponentially for large x and behaves as x~'^'^^ for small x. The hierarchy 
of the correlation functions reasonably agrees with numerical results for the CDM model 
(Bouchet, Schaeffer, & Davis 1992). 

Figure 15 shows the density distribution function multiplied by density {pdn/dp) for 
the 50 Mpc box C+HDM simulation at different redshifts. The distribution of the total 
density (both "cold" and "hot" components included) was obtained on 256^ mesh using 
the Cloud-In-Cell algorithm. We do not show the function at small densities because 
discreteness effects are too large. Figure 15 presents values of Y^- Pi{p'^ p + Ap)/ApN, N = 
256^, which is an estimator for pP{p). The full curves are labeled by the expansion 
parameter a. The function P{p) grows with time and its slope becomes smaller. By 
redshift z = 1 the function has a power-law shape with the slope pP{p) oc p~^, which 
is shown as the dashed line. This is consistent with a large fraction of matter being in 
pancake-type structures. The total density contrast (as estimated by the linear theory) at 
that moment was (5p/p)totai = 1-17, which is the amplitude characteristic for the formation 
of filaments. 

At the final moment 2; = the density distribution showed the power law behavior at 
small densities with the slope ^ —2.3 and an exponential fall off at large densities. The 
dot-dashed curve in Figure 15 shows the fit to the distribution function 

P(p) = 0.275 (A)"'"exp(-^). p>10(p>. (12) 
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Let us compare with the results for the CDM model: Bouchet, Schaeffer, & Davis (1992) 
present counts of particles in cubic cells for a CDM simulation normalized with biasing 
parameter b = 2.5. The size of their computational box was 64 Mpc. If we take the CDM 
results for the size of the cubic cell having on average the same number of particles (1/8) 
as a cell in our model, then the slopes of the distribution function are the same (—2.3) 
in both models. Bouchet, Schaeffer, & Davis also give an approximation for the function 
h{N/Nc) introduced above: h{x) ~ 0.075a:~^ exp(— O.OSs). For our C+HDM simulation 
Nc = ^ ^ 105, = p/(p), and h{x) ^ 0.086x~^-^^ exp(— 0.046a;), which is not much 
different from the results for the CDM simulation. 

9. Conclusions and Discussion 

Let us begin by summarizing the conclusions we have reached above from our C+HDM 
high-resolution PM simulations in 14, 50, and 200 Mpc boxes (always with Hq = 50 km 

Mpc~^), our CDM comparison simulation in a 50 Mpc box, and our Press-Schechter 
calculations. 

1) The correlation function of the dark matter evolves as a power law, with slope increasing 
with time. At redshift z = 0, the slope is 7 = 1.8 and the correlation length is 3h~^ Mpc . 
The correlation length for "galaxies" in the model is (5 — 6)h~^ Mpc and it is almost 
independent of the parameters of the biasing scheme. The "galaxy" correlation function in 
redshift space matches that of the CfA slices (Vogeley et al. 1992) reasonably accurately. 

2) Pair-wise velocities of galaxies in the C+HDM model are in a good agreement with the 
results for the CfA catalog (Davis & Peebles 1983): projected velocity dispersion slightly 
increases with separation and its amplitude at the distance 1 Mpc is (340 — 380) km 
This agreement must be regarded as an important success of the model. 

3) At ^ = the mass distribution of dark halos is roughly approximated by dn/d\ogM oc 
M~^, which is close to what was found for the CDM model, but the amplitude of this 
approximation is about a factor of two lower for the C+HDM as compared to the CDM 
model with the same biasing parameter. The Press- Scliccter approximation (Equations 
7, 8) provides a reasonably good fit for the number of dark halos at all redshifts. At 
higher redshifts the C+HDM model predicts a progressively smaller number of dark halos 
of galactic mass. By redshift z — 2 the difference is about a factor of ten. We note that 
the C+HDM model predicts roughly the same number density of galactic halos as the 
CDM model with the biasing parameter 6 = 2.5. Although the C+HDM model does not 
contradict the number of galaxies and quasars at 2; < (2 — 3), more observational data will 
be needed to test the theory at higher redshifts. 

4) The density distribution function at redshift z = behaves as a power law with an 
exponential falloff. The distribution agrees well with that for the CDM model and with 
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the predictions for the hierarchy of AT— body correlation functions of Bahan & Schaeffer 
(1988). 

It was aheady known from the earher linear C+HDM calculations listed in §1 that 
C+HDM provides a good fit to the available large scale structure data such as galaxy 
large-angle correlations and rich galaxy cluster spatial correlations. Figure 16 shows the 
remarkable agreement between the C+HDM predicted bulk flow velocity as a function of 
radius and the POTENT analysis of all the peculiar velocity data now available (Dekel 
1992). That the same model can also account for the few-Mpc-scale velocity data is a 
major triumph for C+HDM, and reflects the special feature of this model: the fact that 
its hot component is less clustered on small scales than the cold and baryonic components. 

Ostriker and Suto (1990) proposed to evaluate models on the basis of a ratio of the 
streaming velocity to the random velocity dispersion of galaxies, which they dubbed the 
"cosmic Mach number" Ai. This ratio has the virtue of being independent of the bias 
factor (i.e. spectrum normalization), to first order. Unfortunately, the available data do 
not yet permit a reliable estimate of M., but for guidance we will use = 2.2 ± 0.5 
for robs = 14.2/j,~^ Mpc (Groth, Juszkiewicz, & Ostriker 1989). Gaussian smoothing 
with rs = 5h~^ Mpc was used to remove small-scale nonlinear velocities. The sample 
radius robs = 14.2/j,~^ Mpc approximately corresponds to the size of a cubical volume 
Lohs = 58 Mpc (we use h = 0.5 as usual, estimate the radius of a top-hat window 
function and equate the volumes in the cubical and spherical windows; in the numerical 
simulations we estimate the Mach number using a cubic real-space window function). For 
comparison, we use linear theory to estimate the rms velocity of matter inside a volume 
of size Lobs = 58 Mpc. To find the noise velocity below 58 Mpc we use the linear theory 
to get the contribution from wavelengths between 58 Mpc and 50 Mpc and we use the 
numerical simulations to estimate the rms noise velocities due to waves between 50 Mpc 
and 16 Mpc, the later corresponding to the small-scale gaussian filter r^ = 5h~^ Mpc. Our 
result for the CDM model Ai = 1.2 is in good agreement with the nonlinear calculation of 
Suto, Cen, & Ostriker (1992). The C+HDM model gives M = 2.2, which perfectly agrees 
with the data quoted. If instead of L = 58 Mpc we choose L = 50 Mpc, the results would 
be slightly higher: M = IA for CDM and M = 2.6 for C+HDM. 

Although the data is uncertain, as we said, C+HDM probably does much better than 
CDM on the Mach number test. Obviously, this test needs to be refined, and a number of 
additional tests are needed, for example: 

• Galaxy velocities on small scales: Analyses of velocities in the local group and other 
groups, the cosmic virial theorem, and infall arguments typically lead to O 0.2 — 0.4. 
Are these results consistent in detail with the predictions of O = 1 C+HDM? 
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• Power spectra and velocities in spheres compared with data: To improve significantly 
on the linear predictions will require large-scale simulations with high enough resolu- 
tion to identify galaxies. 

• Galaxy, cluster, and quasar formation and number densities, ionization of IGM: This 
could be addressed by high-resolution PM simulations that stop at a redshift of order 
1-2, but eventually will require P^M or tree codes including hydrodynamics and star 
formation. These calculations should also clarify the C+HDM predictions for the 
evolution of the number densities of galaxies, groups, and clusters as functions of 
their masses or characteristic velocities, both in high- and low-density regions. New 
X-ray data on clusters (e.g. Edge et al. 1990, Henry & Arnaud 1991) could allow a 
critical test of the model (see e.g. Kaiser 1991). 

• CMB anisotropics: Current small-angle experiments should be close to detecting tem- 
perature anisotropics if the C-I-HDM model is correct. The compatibility of C-I-HDM 
with published microwave background anisotropy results is discussed by Wright et al. 
(1992). C+HDM predicts anisotropies for the OVRO experiment (Readhead et al. 
1989) and the MAX experiment (Devlin et al. 1992) that are near the upper limits 
derived from these experiments. The model is incompatible with the result of Gaier 
et al. (1992), who measure an upper limit to C(0) of 1.4e-5 from their South Pole 
experiment. We note, however, that this result was derived assuming a Gaussian form 
for the radiation correlation function, which is not a particularly good approximation 
to the function derived for our model; further calculations are required to determine 
whether this makes a significant difference in the derived result. Moreover, this up- 
per limit was obtained considering only the channel with the least variation; a more 
conservative treatment leads to limits that are compatible with C-I-HDM (Dodelson 
&; Jubas 1993; J. R. Bond, private communication). We also note that some fine 
tuning of the model predictions for C-I-HDM can be done by changing the amount of 
baryons in the model; for example, at a scale of 1.2 degrees the predictions for AT/T 
are reduced by about 10% if 0,^ is reduced from 0.10 to 0.05. 

All other potentially viable cosmological models should of course also be subjected 
to detailed confrontation with data. Similar tests should therefore be applied to all the 
main competetors to C-I-HDM, such as f2 = 0.2 CDM, f2 = 1 CDM with designer spectra, 
or = 1 cosmic strings or other non-Gaussian fluctuation models with hot dark matter. 
Based on the results that have been obtained thus far for C-I-HDM, we have reason to 
expect that it will do very well. 
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Figure Captions 

Figure 1— The spectra of fluctuations at different redshifts in the "hot" (top panel) and 
"cold" (bottom panel) components. Curves from top to the bottom correspond to redshifts 
2; = 0, 5, 10, 15, 20, 25. The dashed curves present the approximations (Equations 1) to the 
spectra. 

Figure 2— Distribution of matter in a rectangular slice of 10 Mpc thickness in the 200 Mpc 
box simulation. The slice was chosen to pass through the highest density concentration in 
the simulation. The density concentration is a rich galaxy cluster in the top right corner of 
the slice, (a) "Cold" particles (8% of the total number), (b) Points for which the density is 
higher than ptotai > 10(/Ototai) (1/6 of the points), (c) The distribution of 1015 "galaxies" 
in the slice. "Galaxies" were first preselected as maxima of density, and then those with 
mass larger than 10^^ Mq are presented. The area of a circle representing a "galaxy" is 
proportional to its mass. 

Figure 3— Distribution of "cold" (left panel) and "hot" (right panel) particles in a slice 
of 220 kpc X 7 Mpc x 7 Mpc in the 14 Mpc simulation. We show 85% of "cold" particles 
(1.3 X 10^ particles) and 16.6% of "hot" particles (1.2 x 10^ particles). Note that while aU 
large condensations in the plot are seen both in "cold" and "hot" particles (for example, 
the large group in the left top corner, a few big "galaxies" in the filament, and the filament 
itself), small halos of "cold" particles do not find counterparts in "hot" particles (for 
example, two "galaxies" in the void at the right middle part). The resolution for this 
simulation is 55 kpc, which is about 4 times smaller than the distance between small tick 
marks. 

Figure 4— Evolution of the number density A'"(> M) of dark halos with mass larger than 
some limit (labeled in units of Mq ). (In this and the next two figures, the vertical axis 
is labeled in units of Mpc~^, with Hq = 50 km s~ Mpc~^.) The three bottom unmarked 
curves (mass limits 3 x lO"*^"*^ Mq and higher) are for the simulation with 50 Mpc box size. 
The marked curves correspond to the number of objects in the 14 Mpc simulation. The 
dot-dashed curves presents results for the Press-Schecter approximation. 

Figure 5— Number density of halos per unit volume and per unit logarithmic mass interval 
at the present epoch. The right full curve presents the number of halos in the 50 Mpc box 
simulation and the left full curve is for the 14 Mpc simulation. The dashed curve presents 
the analytical approximation Equation (11). Results for the CDM model are shown as the 
dot-dashed curve. 

Figure 6 — Comparison of the number densities of "galaxies" in the C+HDM and CDM 
models. The triangles represent results of the CDM numerical simulations (box size 
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50 Mpc) and the dashed curves are for the Press-Schecter approximation with the param- 
eter 5c = 1.60. Full curves are the Press-Schecter approximation for the C+HDM model. 
Both models have the same biasing parameter b = 1.5. The dot-dashed curve represents 
the Press-Schecter approximation for CDM with b = 2.5, Sc = 1.60, and M = 10^^ M© . 

Figure 7— Evolution of the correlation function of the total density in the 50 Mpc box 
simulation. The curves in the plot are labeled by the expansion parameter a. At the final 
moment the slope of the correlation function is about 7 = 1.8, but the correlation length 
2.2h~^ Mpc is too small as compared to the correlation length of galaxies. The dashed 
curve in this figure shows the correlation function of the "cold" matter. Results for our 
CDM simulation are shown as the dot-dashed curve. 

Figure 8— Correlation function of "cold" particles at the final moment for the 50 Mpc 
box (bottom dashed curve) and for the 200 Mpc box (bottom full curve) . The top dashed 
and full curves represent correlation functions of "galaxies" for 50 Mpc and 200 Mpc 
simulations correspondingly. The dot-dashed line is the power law ^ — (5.5/i~^/r)^'^. 

Figure 9— Influence of different effects on the correlation function, (a) Simulation with 
50 Mpc box. Full curves are for thresholded density (10 for bottom and 25 for top curves). 
The dashed and dot-dashed curves show results for "galaxies" found as maxima of the 
number of "cold" particles inside a sphere of 100 kpc radius. The dashed curve is for 
"galaxies" more massive than 3 x 10^° M© (913 objects). The dot-dashed curve is for the 
87 "galaxies" with mass larger than 3 x 10^^ Mq . The markers represent the correlation 
function of galaxies in the CfA catalog (Davis and Peebles 1983). (b) Effects of peculiar 
motions on the correlation function. Bottom curves are for the 50 Mpc and top curves are 
for the 200 Mpc simulation. Correlation functions of 913 and 2831 "galaxies" in redshift 
space (real space) are shown as solid (dashed) curves. The markers represent the correlation 
function in redshift space for galaxies in the CfA slices (Vogeley et al. 1992). 

Figure 10— Components of the pair- wise velocities in the C+HDM 50 Mpc box simulation. 
The top panel shows (i'2i) for "cold" particles (dashed curve) and for "galaxies" (full curve). 
The dot-dashed curve shows where the velocity is equal to the Hubble velocity: a pair on 
the curve would have zero proper velocity. The bottom panel presents (7\\ = v\\^rms (full 
curves) and a_i_ = v_i_^rms (dashed curves) for "cold" particles (top two curves) and for 
"galaxies" (bottom two curves). 

Figure 11— Same as Figure 10, but for the CDM model. 

Figure 12— Velocity dispersion of "galaxies" (full curves) along the line of sight as a 
function of projected distance in the 50 Mpc box simulations for the C-f-HDM (bottom 
panel) and CDM (top panel) models. The dashed curves are for all "cold" particles and 
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the dot-dashed curves represent the observed behavior 340(/j.r)°-^^ km s ^. Only pairs of 
objects with AV < 1000 km are taken into account. 

Figure 13 — Dependance of the projected pair-wise velocities on different effects. The 
top panel is for the 200 Mpc box simulation. The solid curve is the velocity dispersion 
of "galaxies" along the line of sight for C-I-HDM, and the dot-dashed curve shows the fit 
a = 380(/j.r)°-^^ km s"-*^, which is compatible with the data. The bottom panel shows results 
for different selection criteria. The long-dashed curve shows results for 481 "galaxies" with 
mass larger than lO"*^-*^ Mq . The short dashed curve presents results for pairs with the 
relative distance smaller than 10h~^ Mpc . The dispersion for the cutoff in the velocity 
difference AV^ < 1500 km s~^ is shown as the dot-dashed curve. 

Figure 14— Distribution of the total velocity of galaxies relative to the rest frame in the 
50 Mpc box simulation (913 objects) as a function of orbital circular velocity. The curve 
presents the averaged rms velocity. The error bars show formally estimated statistical 
uncertainties. 

Figure 15— The density distribution function multiplied by density for the C-I-HDM 
50 Mpc box simulation. Values of pP{p) = J2i Pi{P'i P + ^p)/ ^P^i ^ = 256^ are shown. 
The full curves are labeled by the expansion parameter a. The dashed line corresponds to 
the power law P{p) oc p~^. The dot-dashed curve is the fit Equation (12). 

Figure 16 — Bulk flow velocity vs. top- hat radius R (for Hq = 50) from linear {b = 1.5) 
calculations for C-I-HDM (solid curve) and CDM (dashed), compared with POTENT data 
with Gaussian smoothing radius 1200 km from Dekel (1992). The data points are not 
independant. The error bars on the data points take into account errors of velocities of 
galaxies and those introduced by POTENT. Typically the errors are about 15%. Because 
of statistical fluctuations (the velocity of one sphere is measured) the theoretical prediction 
looks like a wide band. We show the band for the 90% confidence level (with this probability 
a mesurement must be above the lower bound or below the upper bound). The middle 
curves correspond to the mean expected values. 

Figures are available by mail from: Prof. Joel Primack, Institute for Particle Physics, Uni- 
versity of California, Santa Cruz, CA 95064 U.S.A. Email (internet): joel@lick.ucsc.edu; 
Fax:(408) 459 3043. 
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